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ABSTRACT 


A  new  method  of  template  matching  is  proposed  using  graphical  templates.  A 
graph  of  landmarks  is  chosen  in  the  tempiate  image.  All  possible  candidates  for  these 
landmarks  are  found  in  the  data  image  using  robust  relational  local  operators.  A 
dynamic  programming  algorithm  on  decomposable  subgraphs  of  the  template  graph 
finds  the  optimal  match  to  a  subset  of  the  candidate  points  in  polynomial  time.  This 
combination  of  local  operators  to  describe  points  of  interest/landmarks  and  a  graph 
to  describe  their  geometric  orientation  in  the  plane,  yields  fast  and  precise  matches 
of  the  model  to  the  data,  with  no  initialization  required.  In  addition  it  provides  a 
generic  tooibox  for  modeling  shape  in  a  variety  of  applications.  This  methodology 
is  applied  in  the  context  of  T2  weighted  MR  images  of  the  brain. 
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§1  Introduction 


The  main  premise  of  this  paper  is  that  recognition  and  registration  can  be 
achieved  through  graphical  models  describing  the  global  planar  arrangement  of  local 
image  features  or  landmarks.  The  local  features  are  defined  through  very  robust 
and  crude  local  operators  describing  local  pixel  intensity  topographies.  They  are 
not  meant  to  identify  any  complex  components  with  great  accuracy.  It  is  the  planar 
configuration  of  these  features,  described  in  terms  of  a  graph,  which  singles  the  true 
positives  from  the  many  false  positives  which  light  up  for  these  operators. 

Model  registration  is  achieved  by  creating  a  graph  (manually  at  this  point)  of 
triangles  describing  the  global  arrangement  of  the  features,  and  then  using  dynamic 
programming  on  decomposable  subgraphs  to  find  the  optimal  match  to  the  large 
collection  of  feature  candidates  which  light  up  in  the  image.  These  methods  are 
applied  to  continuously  varying  objects  appearing  in  medical  imaging.  In  this  paper 
we  present  the  specific  case  of  MRI  brain  scans.  In  Amit  and  Kong  (1993)  the 
method  was  applied  to  hand  xrays. 

Using  similar  ingredients  in  Amit  and  Geman  (1994)  we  describe  an  object 
recognition  algorithm,  implemented  in  the  context  of  handwritten  character  recog¬ 
nition.  Decision  trees  are  grown,  where  the  questions  at  the  nodes  involve  incre¬ 
menting  graphs  of  binary  relations  between  landmarks.  The  binary  relations  involve 
the  relative  orientation  of  the  two  landmarks.  At  each  node  in  the  tree  more  com¬ 
plex  and  informative  questions  are  constructed  using  questions  from  previous  nodes 
which  were  answered  in  the  affirmative.  The  terminal  nodes  of  these  trees  may  in 
the  future  provide  graphical  shape  models  for  the  corresponding  classes. 

The  problem  of  model  registration  in  medical  imaging  is  becoming  of  growing 
interest.  First  it  provides  an  automatic  means  of  identifying  the  various  components 
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of  the  object  in  the  image  (figure  5).  and  in  some  cases  of  segmenting  the  image. 
Secondly  it  provides  a  means  of  studying  the  variability  in  the  family,  classifying 
subgroups  and  identifying  abnormalities.  Thirdly  it  can  provide  a  means  for  coding 
or  data  compression.  It  is  also  related  to  issues  of  data  fusion. 

The  most  extensively  studied  models  come  under  the  general  title  of  elastic 
matching. 


Elastic  matching  versus  landmark  matching 

Grenander  (1970)  introduced  the  idea  of  elastic  deformations  of  one  and  two 
dimensional  templates.  These  ideas  were  implemented  in  a  Bayesian  framework 
using  spectral  representations  of  the  planar  maps  in  Amit  et.  al.  (1991)  for  hand 
xrays  and  in  Miller  et.  al.  (1993)  for  MRI  images  of  the  brain.  Similar  ideas  were 
proposed  by  Bajcsy  and  Kovacic  (1989)  for  MRI  brain  images,  using  techniques 
developed  in  the  optical  flow  and  image  sequence  analysis  literature.  Horn  and 
Schunck  (1980),  Huang  and  Tsai  (1981)  and  Nagel  (1983),  Terzopoulos  (1988). 
A  comparison  of  these  methods  both  from  the  point  of  view  of  the  optimization 
problem  being  posed  and  the  numerical  techniques  can  be  found  in  Amit  (1994). 
One  dimensional  elastic  models  have  also  been  extensively  studied  see  Grenander 
et.  al  (1991).  and  Kass  et.  al.  (1987). 

There  are  several  limitations  to  these  elastic  models.  First  the  matching  cri¬ 
terion  seeks  to  minimize  the  mean  square  over  all  pixels  of  the  difference  between 
the  intensity  of  the  deformed  template  and  that  of  the  data  image.  On  one  hand 
this  is  a  very  well  defined  criterion  however  it  does  not  ensure  that  specific  points 
of  interest  or  landmarks  be  matched  with  great  precision.  This  later  criterion  is 
very  precisely  defined  for  the  human  eye  but  very  hard  to  formulate  in  mathemat- 
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ical  terms,  second  the  deformations  being  used  are  generic  or  non-parametric  in 
nature  and  do  not  depend  on  the  specific  family  of  objects.  Third,  because  of  the 
inherent  non-linearity  of  the  problem,  and  the  fact  that  the  deformations  are  high 
dimensional,  the  computational  tools  for  calculating  the  match  must  use  relaxation 
techniques  which  run  the  risk  of  converging  to  a  local  minimum  which  corresponds 
to  a  poor  match.  This  is  a  serious  problem  in  the  one  dimensional  elastic  algorithms. 

These  issues  point  to  three  categories  which  must  be  addressed  simultaneously 
and  in  an  integrated  manner  in  any  approach  to  the  template  matching  problem. 
The  data  term  which  drives  the  matching.  The  model  and  its  variability.  The 
computational  tools  and  their  limitations. 

The  principal  tool  in  the  method  described  here  is  a  graphical  model  of  land¬ 
mark  points  which  describes  their  planar  arrangement,  together  with  local  operators 
which  identify  candidates  for  the  various  landmarks  in  the  data  image.  This  is  es¬ 
sentially  a  new  approach  to  the  general  program  outlined  in  Haralick  and  Shapiro 
(1993,  16.1.2).  Namely  image  features  are  defined  and  then  correspondence  is  ob¬ 
tained  through  a  consistency  model ,  which  is  optimized  through  an  efficient  al¬ 
gorithm  consisting  of  dynamic  programming  on  decomposable  subgraphs  of  the 
original  model. 

Thus  more  emphasis  is  put  on  the  precise  matching  of  landmarks,  see  Bookstein 
(1991),  local  features  or  points  of  interest  (see  Haralick  and  Shapiro  (1993)  ).  These 
landmarks  are  both  important  for  understanding  and  analyzing  the  image  and  can 
be  identified  using  various  local  operators  which  employ  more  information  than 
individual  pixel  intensity.  Moreover  modeling  the  variability  in  terms  of  the  relative 
locations  o:  the  landmarks  yields  a  more  specific  and  lower  dimensional  description 
of  the  variability  within  a  certain  family  of  objects. 


o 


From  graph  matching  to  planar  maps:  Once  the  graph  has  been  matched  to 
the  data  it  is  possible  to  create  a  planar  map  through  interpolation  as  in  Bookstein 
'1991).  Furthermore  this  interpolated  map  can  serve  as  an  initial  point  for  the 
elastic  matching  algorithms  where  the  correct  match  of  the  landmarks  is  ensured. 

From  graph  matching  to  component  identification:  Given  the  location  of  vari¬ 
ous  components  of  interest  relative  to  the  landmarks,  these  can  be  automatically 
identified  once  the  graph  has  been  matched.  For  example  the  various  anatomies 
present  in  the  T2  weighted  MRI  brain  scan  in  figure  5  were  identified  automatically 
using  the  match  displayed  in  figure  4. 

Ultimately,  for  each  specific  anatomy  or  image  family  we  will  have  a  graph 
which  will  be  automatically  matched  to  any  incoming  data  image,  yielding  an  au¬ 
tomatic  identification  of  the  various  components.  This  method  can  be  thought  of 
as  a  generic  toolbox  for  modeling  shape. 

In  section  2  some  of  the  details  of  the  graph  matching  algorithm  are  described 
together  with  some  of  the  experiments. 


§2  The  graph  matching  algorithm 

A  template  image  is  chosen  from  the  family  for  the  construction  of  the  local 
operators  and  graphical  model.  For  the  purpose  of  obt  aining  planar  maps  and  image 
compression  this  template  image  is  considered  as  part  of  the  model,  together  with 
the  graph  of  landmarks. 

2.1  The  local  operators 

A  collection  of  landmarks  is  chosen  in  the  template  image  (figure  1 )  and  a 
local  operator  is  chosen  for  each  landmark.  The  local  operators  are  designed  to  be 
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robust,  and  hence  crude,  descriptions  of  the  local  topography  of  the  pixel  intensity 
function  in  the  neighborhood  of  the  landmark. 

A  certain  size  neighborhood  is  chosen,  say  mxn,  together  with  an  array  L  of  1- 
s.  —  1-s  and  0-s  of  those  dimensions.  The  sign  of  the  difference  between  the  intensity 
at  a  given  pixel  i,j,  and  the  intensity  at  each  of  the  pixels  in  its  m  x  n  neighborhood 
is  calculated,  to  yield  an  array  of  1-s  and  —  1-s,  at  pixel  i,j.  If  the  percentage 
of  matched  1-s  and  —  1-s  between  L  and  is  above  a  prescribed  threshold  the 

pixel  is  considered  a  candidate  for  the  corresponding  landmark.  The  0-s  region  in  L 
represents  a  region  where  the  results  of  the  A^1'^  are  ignored.  The  simplest  example 
is  L i  ail  1-s.  and  L-\  all  0-s.  This  corresponds  to  a  local  maximum  operator.  The 
reverse  corresponds  to  a  local  minimum  operator.  The  pixels  corresponding  to  a 
•certain  local  operator  typically  occur  in  clusters.  These  were  identified  and  the 
average  location  calculated. 

In  the  experiments  described  below  on  MRI  brain  scans  the  operators  for  the 
eight  landmarks  were  obtained  by  pointing  to  the  eight  points  in  the  template 
image,  extracting  an  11  x  11  neighborhood  and  calculating  the  matrix  of  signs  of 
differences  with  respect  to  the  middle  pixel.  Below  is  an  example  of  the  operator 
used  to  identify  candidates  for  landmark  no.  1  in  figure  1. 
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The  left  panel  in  figure  2  displays  as  circles  the  candidates  for  this  operator 
using  a  threshold  of  .75  The  right  panel  displays  the  outcome  with  a  threshold 
of  .8.  False  negatives  are  to  be  avoided  at  all  costs  because  then  there  will  not 
be  a  match  in  the  image  for  a  specific  landmark.  This  is  done  by  lowering  the 
thresholds  and  increasing  the  number  of  false  positives.  Dealing  with  these  is  the 
role  of  the  graphical  model,  which  provides  constraints  on  the  planar  arrangement 
of  the  landmarks.  Of  course  too  many  false  positives  will  lead  to  erroneous  matches. 

2.2  Decomposable  graphs  for  matching 

In  figure  1  the  eight  landmarks  chosen  in  an  MRI  scan  of  the  brain  are  displayed. 


The  landmarks  and  the  graph  in  the  template  image 
FIGURE  1. 


T\  picallv  the  same  local  operator  will  be  used  for  several  landmarks,  and  many 
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more  candidates  will  be  found  in  the  data  image,  than  landmarks  associated  with 
that  operator  in  the  model,  see  figure  2. 


Candidates  for  the  first  4  landmarks  found  in  the  data.  1  —  o,  2  —  +.3  —  *,4  —  x. 
Left  panel:  threshold  of  .75.  Right  panel:  threshold  of  .8 


FIGURE  2. 

To  find  the  correct  match  of  the  landmarks  in  the  model  to  the  candidate 
pixels  in  the  image  it  is  necessary  to  introduce  constraints  on  the  relative  locations 
of  the  landmarks  in  the  plane.  This  can  be  done  by  defining  a  collection  of  triangles 
between  triples  of  landmarks  identified  on  the  template  image.  At  present  this  is 
done  manually  by  the  user.  See  figure  1.  A  cost  function  is  associated  with  each 
such  triangie  which  penalizes  its  shape  deviation  from  the  corresponding  triangle  in 
the  template  image.  The  total  cost  function  is  the  sum  of  the  cost  functions  over 
all  the  triangles  in  the  model.  The  collection  of  triangles  can  also  be  expressed  as 
a  colored  graph,  where  the  nodes  are  the  landmarks,  the  color  or  type  of  a  node  is 
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given  by  the  type  of  local  operator  used  to  find  its  candidates.  The  edges  in  the 
graph  exist  between  each  two  landmarks  which  belong  to  a  triangle.  This  graph  is 
called  the  template  graph.  In  this  context  the  deformable  template  is  not  an  image 
but  a  graph  model.  It  is  important  to  note  that  these  constraints  need  not  be  local. 
indeed  they  could  involve  far  apart  landmarks.  This  is  in  contrast  to  the  implicit 
constraints  used  in  the  elasticity  models  which  essentially  penalize  large  changes  in 
the  local  lattice  elements. 

Finding  the  optimal  match  then  reduces  to  an  inexact  consistent-labeling  prob¬ 
lem,  (see  Haralick  and  Shapiro  1 1993.  Chapter  IT),  which  is  genericallv  exponential 
in  complexity.  However  if  the  template  graph  is  chosen  so  as  to  be  decomposable, 
it  is  possible  to  find  the  optimal  match  in  polynomial  time  using  dynamic  program¬ 
ming  on  the  graph,  see  for  example  Rose  et.  al.  (1976).  Decomposability  in  the 
present  context  means  that  there  exists  an  order  in  which  the  triangles  of  the  graph 
can  be  successively  eliminated,  such  that  each  triangle  in  its  turn  has  a  free  vertex 
contained  in  no  other  triangle.  W  hen  the  free  vertex  and  the  two  edges  emanating 

from  it  are  removed,  one  of  the  vertices  of  the  next  triangle  in  the  order  is  freed 
and  so  on. 


Dynamic  programming  on  the  graph 

The  graph  displayed  in  figure  1  is  decomposable,  and  the  order  of  elimination 
coincides  with  the  index  number  of  the  landmark  type.  Let  v,,  i  =  1 . 8  de¬ 

note  the  vertices  in  the  graph.  Z,  the  corresponding  local  operator,  and  da  denote 
candidate  points  found  in  a  data  image.  For  each  triangle  i.j,k  in  the  graph  let 
<t>ijk(dQi,dQ; .  dak )  denote  the  cost  incurred  by  matching  v,  to  da{,  vj  to  dQ.  and  vk 
to  dak.  Obviously  dQ{  was  identified  with  the  operator  L,  etc.  The  cost  of  a  full 
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match  is  simply 


^  ^  Qijkidoti  ,  daj ,  dak ). 

(i,j,k)ec 

where  C  is  the  collection  of  triangles  in  the  graph. 

Dynamic  programming  in  this  context  is  implemented  as  follows.  For  each 
possible  match  da3,da 3  for  vertices  V2.V3  find  the  point  dax  which  minimizes 

©123  (dQ  1 , da3 ,  da3 ), 

call  the  minimizing  index  01(0:2  ,<*3)  and  the  associated  cost  £123(02,03).  Store 
these  two  for  each  pair  (02,03).  Now  for  each  possible  match  daa,daA  find  the 
point  da,  which  minimizes 

<?234 (da-2  1  da3,  da4  )  +  £123(02,03). 

Note  that  the  second  term  is  already  stored  in  memory.  Again  call  the  minimizing 
point  02' 03.04)  and  the  associated  cost  £234(03,04),  and  so  on  until  the  best  dae 
is  found  for  every  pair  of  daT,dQa  with  the  updated  cost  function 

£678  (  O7 .  Og  )  0(37g  (  d(^3 ,  d^j ,  d o-g  )  £567 (og  ,07). 

The  optimal  pair  07,  o§  is  chosen,  Og  is  taken  to  be  the  already  stored  06(07, 03),  for 
05  take  o 5 ( o g ,  07)  and  so  on,  until  the  entire  optimal  match  aj, . . . ,  o|  is  recovered. 

In  abstract  terms  the  graph  shown  in  figure  1  is  simply  a  second  order  Markov 
graph.  This  is  not  always  the  case.  In  Amit  and  Kong  (1993)  more  complex  graphs 
were  required.  Even  in  the  present  context  replacing  the  edge  between  vertices  3 
and  5  to  an  edge  between  vertices  3  and  6  yields  comparable  results,  however  the 
graph  is  no  longer  second  order  Markov.  The  match  obtained  for  the  candidate 
points  shown  in  the  left  panel  of  figure  2  is  shown  in  figure  3.  For  vertex  8,  the 
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The  match  to  the  candidates  in  the  left  panel  of  figure  2 


FIGURE  3. 

algorithm  prefered  a  candidate  slightly  to  the  left  of  the  correct  one.  so  not  all  is 
perfect.  We  address  this  issue  in  the  discussion  section.  The  match  obtained  for 
the  candidate  points  shown  on  the  right  panel  of  figure  2  is  shown  in  the  left  panel 
of  figure  4  together  with  the  match  for  another  brain  scan  displayed  on  the  right. 

In  figure  5  the  automatic  identification  of  anatomies  obtained  from  the  graph 
match  is  displayed. 

Dynamic  programming  has  been  used  in  image  analysis  in  a  variety  of  contexts 
such  as  road  tracking  (Barzohar  and  Cooper  (1993))  ,  stereo  (Ohta  and  Kanade 
(1985)),  and  artery  tracking  (Petrocelli  et.  al.  (1992)).  All  these  settings  are  one 
dimensional  in  nature,  and  the  constraints  enforced  by  the  underlying  graph  are 
all  local.  To  our  knowledge  this  is  the  first  attempt  to  incorporate  the  efficient 
computational  tool  of  dynamic  programming  in  an  inherently  two  dimensional  and 
non-local  imaging  problem. 
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Matches  for  two  different  images,  with  candidates  obtained  with  the  .8  threshhold. 


FIGURE  4. 

Decomposabilitv  imposes  certain  limitations  on  the  graph,  and  may  limit  the 
geometric  and  topological  information  it  contains.  In  Amit  and  Kong  (1993)  this  is 
addressed  by  starting  from  a  non-decomposable  graph  which  contains  the  desired 
information,  splitting  it  into  decomposable  subgraphs  with  certain  points  in  com¬ 
mon.  An  optimal  match  is  found  using  the  first  graph.  When  the  optimization 
procedure  begins  for  the  second  subgraph  the  match  for  this  subset  of  data  points 
is  fixed  as  the  match  provided  by  the  first  subgraph.  Now  the  optimal  match  is 
found  for  the  remaining  vertices  of  the  second  subgraph.  In  the  problem  presented 
here  one  decomposable  graph  was  sufficient  to  obtain  the  matches. 
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Automatic  identification  of  anatomies  in  the  two  images. 


FIGURE  5. 


The  cost  functions 

The  cost  function  o  has  the  form 


4*ijk  (  i  dQ ,■  ,  dQk  ) 


l°g  (K^ati i  dQk  )/l(da; ,  dQ. )) 
log  (l(daj ,  dQk  )/l(da. ,  dai)) 


+ 


log  (l{da„dai)/l(dak,da.)) 


log  (l(vi,Vk)/l(Vi,Vj)) 
log  {l(Vj,Vk)/l(vj,Vi)) 
log(l(vk,Vj)/l(vk,  V,)) 


where  /  denotes  the  Euclidean  distance  between  between  two  points.  In  other  words 
the  sum  of  the  square  of  the  difference  of  the  log  ratios  of  the  lengths  of  each 
pair  of  edges  between  the  candidate  triangle  da. ,  da. ,  dai  and  the  template  triangle 

v*'vi'  l'k-  In  addition  a  hard  constraint  was  imposed  which  did  not  allow  an  angle 
to  change  by  more  than  7t/4. 

Note  that  these  cost  functions  are  completely  scale  and  rotation  invariant.  On 
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the  other  hand  the  local  operators  although  robust  to  large  changes  in  scale  are 
robust  to  small  changes  in  rotation.  In  figure  6  we  show  a  successful  match  to  a 
rotated  version  (7 r/8)  of  the  brain  scan.  Similar  results  were  found  for  the  other 
images,  within  a  range  of  (-tt/S. 7t/8). 


A  match  obtained  for  an  image  rotated  by  7r/8. 


FIGURE  6. 

At  tt/4  the  algorithm  failed,  the  operators  were  completely  missing  the  correct 
points.  Note  that  a  match  of  the  graph  provides  the  orientation  of  the  object  as 
a  side  benefit.  This  partial  rotation  invariance  for  a  certain  range  of  angles  in  the 
neighborhood  of  zero,  can  yield  a  fully  rotation  invariant  algorithm.  In  the  present 
context  the  image  would  be  rotated  at  0, tt/4, 7t/2, ...  and  for  each  such  angle  a 
match  would  be  obtained.  The  angle  with  a  match  of  minimal  cost  would  be  used 
to  identify  the  correct  match  and  the  correct  orientation. 
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Computation 


The  landmarks  candidates  are  identified  in  the  data  through  some  routines 
operated  through  MATLAB.  and  take  10-15  seconds  on  a  SUN  SPARC  10.  The 
software  for  obtaining  the  optimal  match  is  written  in  C.  and  also  takes  10-15 
seconds  on  the  same  machine  with  approximately  400  candidates  all  together  for 
all  landmarks,  (about  50  per  landmark).  The  software  is  available  upon  request. 


§3  Discussion 

The  strength  of  this  algorithm  is  its  computational  efficiency  and  the  absence 
of  any  need  to  initialize  the  matching  algorithm  or  the  optimization  procedure.  This 
is  one  of  the  drawbacks  of  many  relaxation  techniques.  If  the  template  or  model  is 
not  initially  set  nearby  the  correct  location  the  relaxation  algorithm  may  end  up  in 
completely  erroneous  solutions.  This  is  for  example  the  problem  with  the  dedicated 
template  models  of  Yuille  et.  al.  (1989),  or  in  Cootes  and  Taylor  (1992b),  as  well 
as  for  the  elastic  matching  models  mentioned  above. 

The  methodology  presented  here  is  new  and  many  question  remain  open.  At 
the  moment  the  masks  for  the  local  operators  are  automatically  obtained  once  the 
landmarks  are  pointed  out  by  the  user  in  a  chosen  template  image.  How  should  the 
thresholds  be  determined,  should  there  be  a  score  associated  with  each  candidate 
according  to  some  criterion?  Such  scores  could  be  estimated  using  training  images 
together  with  some  statistic  of  the  pixel  configuration  in  the  neighborhood  of  the 
landmark.  This  may  help  resolve  local  mismatches  as  in  figure  3.  How  would 
other  types  of  the  many  local  operators  proposed  in  the  computer  vision  literature 
perform  in  conjunction  with  the  graph  models?  Could  different  operators  be  used 
at  different  image  resolutions  to  obtain  matching  at  various  levels  of  detail? 
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There  remains  the  interesting  question  of  how  to  automate  the  creation  of  the 
graph.  It  is  hard  to  imagine  a  fully  automated  mechanism,  because  in  many  cases 
decomposaole  graphs  will  be  insufficient.  It  is  then  necessary  to  iterate  between  de- 
composaDie  subgraphs.  However  automatic  generation  of  trial  decomposable  graphs 
with  certain  constraints  or  conditions  should  be  possible  to  implement,  thus  aiding 
the  user  in  developing  the  full  graph  to  be  used  in  the  algorithm. 

Finally  the  present  optimization  algorithm  only  identifies  the  best  match.  Find¬ 
ing  the  best  k  matches  should  make  the  procedure  more  robust,  and  is  still  feasible 
in  polynomial  time.  An  interesting  problem  would  be  how  to  chose  from  among 
these  k  best  matches. 

Another  interesting  issue  is  which  image  to  use  to  generate  the  model.  Can 
the  model  be  improved  as  new  images  come  in  and  are  matched?  For  example 
certain  parameters  of  the  cost  function  could  be  updated  to  reflect  the  individual 
distributions  of  the  triangles  in  the  graph.  Even  the  ‘mean’  triangle  provided  by 
the  initial  graph  model  may  change. 
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